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Abstract.  We  present  methodology  for  obtaining  forward  solutions  to  Maxwell’s 
equations  in  two  dimensions,  in  the  presence  of  a  Debye  medium.  Perfectly  Matched 
Layer  (PML)  absorbing  boundary  conditions  are  used  to  absorb  incoming  energy  at  the 
finite  boundaries.  A  time-domain,  PDE  formulation  is  presented,  and  a  finite  difference 
time-domain  (FDTD)  algorithm  is  used  to  obtain  numerical  solutions.  A  least  squares 
formulation  of  the  inverse  problem  results  from  a  careful  consideration  of  the  noise 
model  for  data  generation.  The  inverse  problem  is  solved  with  varying  levels  of  noise 
in  the  data,  and  a  frequency  domain  analysis  is  given  that  provides  an  explanation  of 
the  results.  The  results  and  analysis  motivate  strategies  for  solving  the  inverse  problem 
that  decrease  computational  cost.  Finally,  a  result  from  the  statistical  theory  of  large 
samples  is  used  to  obtain  estimates  of  the  variability  in  parameter  estimates  that  is 
due  to  the  variability  in  the  noisy  data. 
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Figure  1.  Infinite  strip  antenna  and  dielectric  slab. 


1.  Introduction 

In  [3]  the  authors  developed  techniques  for  solving  Maxwell’s  equations  in  one  dimension 
in  a  Debye  medium,  and  for  solving  the  associated  parameter  estimation  or  inverse 
problem.  The  experimental  setup  assumed  in  the  formulation  of  the  one  dimensional 
problem  includes  an  antenna  that  has  infinite  extension  in  both  the  x  and  y  directions 
(i.e.,  an  infinite  sheet).  In  this  paper,  we  consider  instead  the  case  of  an  infinite  strip 
antenna  which  is  infinite  in  the  y  direction,  but  is  finite  in  the  x  direction.  This  setup 
is  illustrated  schematically  in  Figure  1.  The  resulting  Maxwell  system  involves  oblique 
incident  waves  on  the  dielectric  slab  and  hence  is  two  dimensional.  The  goal  of  this 
paper  is  to  extend  results  in  [3]  to  this  problem. 

Regardless  of  the  experimental  setup,  in  practice  absorbing  boundary  conditions 
are  necessary  in  order  to  model  the  (effectively)  infinite  spatial  domain  by  a  finite 
computational  domain.  In  the  one  dimensional  case,  perfectly  absorbing  boundary 
conditions  exist  and  are  easily  implemented  [3].  Unfortunately,  such  boundary 
conditions  do  not  exist  in  higher  dimensions.  Though  various  analytic  ’’absorbing” 
boundary  conditions  are  available  for  use  on  the  two  dimensional  problem,  a  more 
effective  technique  employs  perfectly  matched  layers  (PMLs)  as  a  fictitious  absorbing 
layer  surrounding  the  region  of  interest  (see  Figure  2).  If  the  PML  is  constructed 
properly,  there  is  no  reflection  at  the  PML  interface  and  any  energy  returning  to  the 
domain  of  interest  after  travelling  through  the  PML  is  negligible.  In  Section  2,  we 
present  a  system  of  PDEs  that  models  the  setup  illustrated  in  Figure  1  with  a  Debye 
dielectric  slab  and  PML  absorbing  boundaries.  This  system  assumes  that  the  electric 
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Figure  2.  Problem  geometry 


field  generated  by  the  antenna  is  polarized  so  that  it  only  has  a  y  component  and  is 
solved  using  a  finite  difference  time-domain  (FDTD)  algorithm  [16].  Numerical  results 
from  forward  simulations  are  presented. 

In  Section  3  we  discuss  the  solution  of  the  associated  parameter  identification 
problem.  That  is,  given  data  generated  from  the  experimental  setup  illustrated  in 
Figure  1  with  a  Debye  dielectric  slab,  determine  the  parameters  that  characterize  the 
Debye  medium.  We  begin  the  section  with  the  noise  model  for  data  generation  for  such 
experiments.  Using  this  noise  model,  we  arc  led  to  the  negative  log- likelihood  function, 
which  takes  a  least  squares  form.  The  inverse  problem  is  then  to  minimize  this  function. 
We  discuss  the  modified  Levenburg-Marquardt  method  that  we  used  to  this  end.  Several 
attempts  are  made  at  the  inverse  problem  with  varying  levels  of  noise.  The  results  are 
explained  and  new  approaches  are  motivated  by  a  frequency  domain  analysis.  This 
analysis  allows  one  to  obtain  a  knowledge  of  which  parameters  are  or  are  not  likely  to 
be  reconstructed  accurately. 
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Because  our  data  is  noisy,  it  can  be  viewed  as  a  realization  of  a  random  variable 
and  is  therefore  variable.  That  is,  in  the  context  of  the  setup  of  Figure  1,  if  we  perform 
exactly  the  same  experiment  many  times,  we  will  obtain  different  data  sets  each  time. 
Consequently,  if  we  solve  the  inverse  problem  with  each  of  these  different  data  sets,  we 
will  obtain  different  values  for  the  reconstructed  parameters.  Thus  the  variability  in 
the  data  translates  into  variability  in  the  parameter  estimates.  In  Section  4,  we  use 
the  statistical  theory  of  large  samples  [9]  to  obtain  estimates  for  the  variability  in  our 
estimates.  Conclusions  are  presented  in  Section  5 


2.  The  Forward  Problem 


We  begin  with  Maxwell’s  equations  for  a  region  with  free  charge  density  p  =  0,  which 
are  given  by 

V  D  =0  (1) 

V  •  B  =  0  (2) 

V  x  E  =  -dtB  (3) 

V  x  H  =  dtD  +  J,  (4) 

where  the  vectors  in  (l)-(4)  are  functions  of  position  r  =  (. x,y,z )  and  time  t,  and 
J  =  Jc  +  Js,  where  Jc  is  the  conduction  current  density  and  Js  is  the  source  current 
density.  We  assume  only  free  space  can  have  a  source  current,  and  thus  either  Jc  =  0 
or  Js  =  0,  depending  on  whether  or  not  the  region  is  free  space.  For  a  Debye  medium, 
magnetic  effects  are  neglected,  and  we  assume  that  Ohm’s  law  governs  the  electric 
conductivity.  Hence,  within  the  dielectric  medium 

B  =  fi0  H,  (5) 


Jr  =  a  E. 


(6) 


The  displacement  vector  Dr  on  the  other  hand,  has  a  more  complex  representation, 
namely, 

D  =  eoeooE  +  P,  (7) 

where  P  is  the  electric  polarization  vector  that  satisfies  the  differential  equation 


tP  +  P  —  eo(o  —  e0C)E 


(8) 


inside  of  the  dielectric  and  is  0  outside  of  the  dielectric. 

The  computational  domain  for  our  problem  is  given  in  Figure  2.  Notice  first  that 
it  is  contained  in  the  x,z  plane.  This  is  facilitated  by  the  fact  that  we  have  assumed 
uniformity  along  the  y- axis.  Secondly,  note  that  the  computational  domain  has  finite 
support.  This  is  desirable  since  in  order  to  numerically  solve  Maxwell’s  equations,  a 
finite  computational  grid  is  necessary.  For  our  problem,  the  resulting  computational 
boundaries  will  generally  reflect  electromagnetic  waves,  and  hence,  numerical  solutions 
will  contain  non-physical  artifacts.  One  of  the  more  effective  ways  to  combat  this 
difficulty  is  to  use  Perfectly  Matched  Layers  (PMLs).  The  PML  is  a  fictitious  medium 
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that  is  placed  on  the  outside  of  the  region  of  interest  (see  Figure  2)  in  order  to  absorb 
incoming  energy  without  reflection. 

In  order  to  be  more  explicit  about  the  geometry  of  the  computational  domain,  we 
introduce  some  notation.  Let  T>  denote  the  computational  domain  as  pictured  in  Figure 
2.  We  partition  the  interval  X  into  disjoint  closed  intervals  X_,  X0,  and  X+  such  that 

max  A"  =  min  X0 
max  X(j  =  min  X+, 

and  partition  the  interval  Z  into  disjoint  closed  intervals  Z_,  Zv,  and  Zd  such  that 
max  Z_  =  min  Zv 
max  Zv  =  min  Zd. 

The  regions  Zv  and  Zd  are  the  vacuum  and  dielectric  regions  respectively.  Let 
Zq  =  Zv  U  Zd- 

The  region  X0  x  Z0  will  be  our  domain  of  interest.  The  buffer  region  T>  \  (X0  x  Z0) 
outside  our  domain  of  interest  contains  the  PMLs.  Note  that  the  PMLs  surround  the 
computational  region  on  only  three  sides.  This  suffices  since  it  is  assumed  that  the 
dielectric  is  backed  by  a  reflective  material  such  as  metal,  where  the  boundary  condition 
E  =  0  is  used. 

Our  source  term  JA.  models  a  finite  antenna  in  free  space.  It  will  have  the  form 

fx(t) 

3s(t,x,z)  =  I[-X,x\{x)  ■  Ss(z)  ■  fy(t )  (9) 

_  fz(t )  _ 

where  P[a,b\  is  the  indicator  function  on  the  interval  [a,  b],  [— x,x]  G  X0,  z  is  the  2- 
coordinate  of  the  position  of  the  antenna.  The  functions  fm,  for  m  =  x,y,z,  are 
determined  by  the  polarization  and  intensity  profile  of  the  interrogating  electromagnetic 
pulse. 

We  now  express  Maxwell’s  equations  in  a  form  that  is  useful  for  our  particular 
problem.  Combining  the  curl  equations  (3)  and  (4)  with  equations  (5),  (6),  (7),  and  (8), 


using  the  rescaling 

E=./^E,  (10) 

V  A*o 

and  writing  the  resulting  equations  in  the  frequency  domain,  we  have 

jivD  =  c0  ■  V  x  H  —  Js  (11) 

D(u;)  =  e*(u)  ■  E(u;)  (12) 

juH  =  —  c0  •  V  x  E,  (13) 


where  denotes  the  Fourier  transform  in  the  time  variable,  and  e*M  is  given  by 


e  (cn)  —  Coo  +  - — | — : - h  - • 

l+jcnr  jue0 


(14) 
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Note  that  the  conductivity  term  Jc  =  <rE  is  now  contained  within  (12)  rather  than 
within  the  forcing  term  in  (11),  as  was  the  case  in  the  original  formulation  (l)-(4). 
In  particular,  it  appears  in  the  last  term  of  (14).  The  polarization  equation  (8)  is 
also  contained  with  (12).  It  is  found  in  the  second  term  of  (14).  We  express  Maxwell’s 
equations  in  this  way  for  three  reasons.  First,  in  this  formulation,  the  parameters  that  we 
hope  to  identify  in  the  inverse  problem  are  contained  entirely  within  the  single  equation 
(12)  with  (14).  Secondly,  given  particular  values  for  the  parameters  e^,  es,  r,  and  a,  for 
a  fixed  frequency  u,  equation  (14)  will  reveal  the  sensitivity  of  the  displacement  vector 
D  to  changes  in  these  parameters.  This  information  will  be  particularly  useful  when 
we  attempt  the  parameter  identification  problem.  Third,  this  formulation  results  in  a 
simplified  analysis,  PML  formulation,  and  PML  implementation. 


2.1.  Reduction  to  Two  Dimensions  and  to  the  TM  Mode 


Since  the  problem  geometry  does  not  depend  on  y,  Maxwell’s  equations  can  be  simplified. 
In  addition,  we  make  the  assumption  that  our  pulse  is  polarized  so  that  it  only  has  a  y 
component.  We  define  fy  in  (9)  by 


fy(t )  =  l[o,tf](t)sin(ut), 


where  c otf  is  an  integral  multiple  of  27t,  and  hence, 


j  s(t,X,z) 


0 

T[-X,x]{x)  ■  Sg(z)  ■  T[o,tf\{t)  ■  sin(ut) 

0 


The  curl  equations  (11), (13)  then  have  the  equivalent  time-domain  formulation 


dtHx(x,z) 

1 

tqi 

_ i 

dtHy(x,z) 

=  c0 

j  dzExiXi 

_  dtHz(x,z)  _ 

&xEy{x ?  Z ) 

dtDx(x,z ) 

- dzHy(x,z ) 

dtDy(x,z) 

=  Co 

dzHx{%',z')  dxHz(x,z) 

i 

v  N 

IQ 

- 1 

dxHy  (x,  z^ 

(15) 


(16) 


(17) 


(18) 


where  Js  is  given  by  (16). 

Assuming  zero  initial  conditions  for  E  and  H.  the  vector  system  (17)-(18)  reduces 
to  the  set  of  three  equations 


dDy  =  (dHz  _  ()H,\ 
dt  C°  \  dx  dz  ) 
dHx  _  dEy 
dt  C°  dz  ’ 
dHz  _  _  dEy 
dt  C°  dx 

which  is  known  as  the  TM  mode. 


J3,y(t,X,  Z) 


(19) 

(20) 
(21) 
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We  choose  to  work  with  the  TM  mode  since  (11)-(13)  then  reduces  to  four  equations. 
If  we  chose  instead  to  polarize  our  electromagnetic  pulse  so  that  JStV  =  0,  while 
Js,x,Js,z  ~f~  0,  the  TE  mode  would  result,  and  (11)-(12)  would  reduce  instead  to  five 
equations.  In  the  TM  case,  (12)  is  needed  only  for  updating  Ey,  whereas  in  the  TE  case 
it  is  used  to  update  both  Ex  and  Ez. 


2.2.  Perfectly  Matched  Layers 


To  simplify  notation  in  the  sequel,  we  suppress  the  notation  found  in  (19)-(21). 
In  order  to  give  physical  motivation  for  the  use  and  effectiveness  of  the  PML  we  first 
consider  (19).  We  note  that  JS)V  =  0  inside  of  the  PML  regions.  To  model  a  lossy 
medium  for  z  6  Z_  (see  Figure  2),  we  add  a  loss  term: 


dDy  =  / dip_dHx\  _  a(z) 

dt  C°  \  dx  dz  J  e0 


(22) 


where  cr(z)  =  0  for  z  G  Z0.  Inside  of  the  PML,  i.e.  for  z  G  Z_,  a(z)  is  a  smooth  function 
of  z  that  increases  from  a  value  of  zero  at  the  free-space/PML  interface  to  a  value  of 
&max  at  the  boundary.  In  the  frequency  domain,  equation  (22)  can  be  written  as 


P  cr[z ) 
ju  ■  ll  +  - - 

V  jue0/ 


Dy  ~  Cq 


'OIL,  dll,. 


dx 


dz 


Making  similar  choices  for  (20)  and  (21),  we  obtain 


V  jueo, 

(  cr{z)' 

V  juto, 

(  d(z)' 

JU)  ■  ll  +  - 

V  Jueo, 


Dy  Cq 


'dip  dHr 


dx  dz 


rr  _  dEy 

Hx  Co  „  , 

dz 

77  _  dEv 

dx 


(23) 


(24) 

(25) 

(26) 


where  the  over  Ey,  Dy,  HX1  and  Hz  denotes  the  Fourier  transform  in  time  variable. 

For  z  E  Z_,  the  PML  is  modelled  in  the  frequency  domain  by  (24)-(26).  These 
very  specific  choices  for  the  fictitious  complex  permittivity  ey(z)  =  1  +  a(z)/juj€o  and 
permeabilities  p>x(z)  =  (1  +  a(z) / joueo)^ ,  p>z(z)  =  1  +  a(z)/ jcue0  result  in  zero  reflection 
(in  the  continuum  space)  of  electromagnetic  energy  at  the  free-space/PML  interface 
[15].  To  see  this,  we  note  that  if  a  wave  is  propagating  in  a  medium  A  and  impinges 
upon  a  medium  B,  the  amount  of  reflection  is  determined  by  the  respective  impedances 
and  is  given  by 


r  =  — — — .  (27) 

V  A  +  rjR 

The  impedance  7/  of  a  medium  with  complex  permittivity  e  and  complex  permeability 
fi  is  given  by 


f)  = 
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For  our  formulation,  e,  /t,  and  rj  are  diagonal  tensors  of  the  form 


e  =diag(ex,ey,ea), 

fl  diag(/Xa;,  fly  ,  fj'z)  1 

V  =  diag  (fjx,f}y,fjz). 


The  dielectric  constants  ey,  fix ,  and  fiz  are  chosen  so  that  they  satisfy  two  sets  of 
conditions,  one  of  which  is  that  the  resulting  diagonal  tensor  T  defined  by  (27)  is  the 
zero  tensor.  For  details,  see  [14,  15]. 

In  discretized  space,  a  perfectly  matched  interface  does  not  exist  [6],  and  hence, 
reflections  do  occur  in  computations.  Fortunately,  they  can  be  made  to  be  negligible. 
In  our  implementation,  this  is  done  by  choosing  the  depth  of  the  PMLs  carefully. 

We  also  need  to  model  the  PMLs  found  in  regions  where  x  G  X_  U  X+.  Using 
arguments  equivalent  to  those  needed  to  derive  (24)-(26),  the  following  system  of 
equations  is  obtained: 


jv 

ju 

ju 


Dy  —  Co 


'dHz 

dx 


77  _  dEv 

Co  „ 
oz 

rr  _  9Ey 

Hz  Co 

ox 


dHx\ 

dz  ) 


(28) 

(29) 

(30) 


This  system  together  with  (12)  gives  us  the  full  system  of  equations  that  we  wish  to 
solve.  (For  a  thorough  overview  of  PML  technology  for  electromagnetics  problems  see 
[16]  and  [8].) 

We  now  must  write  (28)-(30)  in  the  time  domain. 


2.3.  The  Full  Formulation  in  the  Time  Domain 

It  is  a  straightforward  exercise  to  show  that  (28)-(30)  can  be  expressed  in  the  time 
domain  as 


dtD*  +  r,!  =  Co  (dxHz  -  dzHx) , 

e0 

dtDy  +  °^~Dy  =  dtD* 

Co 

d,Hl  +  —h:  =  Co  d,E„ 

Co 

dtHx  =  dtH*x  + 

Co 

dtH*  +  ^H*z  =  -  c0dxEy, 

Co 

dtHz  =  dtHl  + 

Co 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 
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where  D*,  H*,  and  H*  are  defined  by  equations  (32),  (34),  and  (36)  respectively.  Note 
that  outside  of  the  PML  regions,  a(x)  =  a(z)  =  0  and  (31)-(36)  reduces  to  (19)-(21). 
(Recall  that  we  have  suppressed  the  notation.) 

The  y  component  of  the  time  domain  expression  of  (12)  is  given  by 

Dy(t )  =  CooEyit)  +  Py(t )  +  Cy(t ),  (37) 


where  the  polarization  Py  satisfies 

Py  T  Py  (U  ^oo)Pyi 
and  the  conductivity  term  Cy  satisfies 
Cy  =  (v/e0)Ey 

inside  of  the  dielectric.  Outside  of  the  dielectric,  Cy  —  Py  —  0. 


(38) 

(39) 

Recall  that  we  rescaled 


Ey  via  equation  (10)  and  then  suppressed  the  corresponding  notation  change.  This 
accounts  for  the  change  in  form  from  (8)  to  (38).  Also,  note  that  Dy  in  (37)  is  not  the  y 
component  of  the  displacement  vector  D  of  equation  (7).  It  is,  rather,  the  y  component 
of  D  in  (12)  written  in  the  time  domain.  The  was  suppressed  at  the  beginning  of 
Section  2.2  in  order  to  simplify  notation. 

Equations  (31)-(38)  are  the  full  set  of  equations  that  we  solve  numerically.  We 
use  finite  differences.  In  particular,  we  use  the  FDTD  algorithm  of  [15,  16],  which  uses 
forward  differences  to  approximate  the  time  and  spacial  derivatives.  The  Yee  staggered 
grid  is  used  [15,  16],  which  ensures  second  order  accuracy  within  the  free  space  region. 
Unfortunately,  this  accuracy  is  not  maintained  within  or  at  the  boundary  of  the  dielectric 
medium.  The  FDTD  algorithm  is  fully  explicit,  and  does  not  require  the  inversion  of 
large,  sparse  matrices.  Our  implementation  of  the  FDTD  algorithm  for  this  problem, 
including  the  implementation  of  the  PMLs,  follows  details  found  in  [15]  closely,  with  the 
exception  that  the  corner  regions  in  which  there  is  overlap  of  PML  layers  are  dealt  with 
as  in  [7].  Gauss’  divergence  laws  (1)  and  (2)  are  automatically  enforced  by  the  FDTD 
algorithm  (see  [16]).  Wcll-posedness  for  this  system  is  an  important  question,  which  we 
hope  to  address  in  a  future  paper. 

In  all  simulations,  the  spatial  and  time  steps  are  related  via 


a  A x 

At=*7 


(40) 


where  co(=  3  x  108)  is  the  (approximate)  speed  of  light.  This  satisfies  the  CFL 
condition,  and  hence,  guarantees  the  stability  of  the  FDTD  algorithm  in  free  space 
[16].  Fortunately,  stability  is  retained  when  a  Debye  medium  is  introduced  [12], 

An  unfortunate  side  effect  of  numerical  solutions  to  Maxwell’s  equations  is  that 
dispersion  is  introduced  by  the  numerical  scheme.  This  is  the  case  for  both  finite 
difference  and  finite  element  (FEM)  methods  (see  [8,  16]).  Techniques  can  be  used  to 
minimize  dispersion  but  such  techniques  are  very  difficult  to  implement  when  dispersive 
materials  are  present  [12],  For  the  Debye  medium  in  particular,  it  is  noted  in  [12]  that 
when  the  FDTD  algorithm  is  used  spurious  dispersion  may  occur  unless  At  ~  (9(1CT3t), 
which  can  be  a  highly  restrictive  condition. 
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A  different,  but  standard,  approach  used  in  solving  Maxwell’s  equations  in  the 
time  domain  is  to  first  construct  the  wave  equation  for  the  electric  held  E,  and  then 
solve  it  numerically.  In  [2]  this  is  done  for  the  forward  problem  addressed  in  this  paper 
using  FEM.  The  main  drawback  of  the  FEM  approach  is  that  at  each  time  step  a  large 
sparse  linear  system  must  be  solved.  This  is  not  the  case  for  the  FDTD  algorithm. 
Consequently,  for  the  same  grid  spacing,  FDTD  is  much  faster.  In  addition,  the 
implementation  of  the  FDTD  algorithm  is  much  simpler  for  this  problem  than  it  is  for 
FEM.  This  is  particularly  true  for  the  implementation  of  the  PMLs.  In  the  numerical 
comparisons  that  we  have  done  between  these  two  approaches,  the  FDTD  algorithm 
was  at  least  ten  times  faster  when  the  same  computational  grid  was  used.  On  the  other 
hand,  we  have  not  carried  out  a  convergence  comparison  between  FDTD  and  FEM  for 
the  forward  problem.  It  is  likely  that  the  forward  solutions  obtained  using  FEM  will 
converge  on  a  coarser  grid  than  does  FDTD,  since  FDTD  loses  its  second  order  accuracy 
within  the  dielectric  and  at  the  diclcctric/free  space  interface.  If  this  is  the  case,  savings 
in  computational  time  may  be  lost.  Hence,  it  is  not  clear  at  this  point  which  approach 
(FEM  or  FDTD)  is  better  for  use  in  solving  the  inverse  problem.  These  questions  are 
currently  being  addressed  in  further  efforts. 

2-4-  The  Forward  Problem  for  a  Debye  Medium 

In  this  section  we  present  details  and  results  from  an  attempt  at  solving  (31)- (38)  using 
the  FDTD  approach  discussed  in  the  previous  section.  We  begin  with  the  geometry 
of  the  computational  region.  The  necessary  measurements  are  given  by  Xq  =  [0,0.1], 
Zv  =  [0,0.2],  Zd  =  [0.2,  0.4]  (these  measurements  are  in  meters).  In  order  to  achieve 
convergence  of  the  finite  difference  approximations  we  take  the  number  of  nodes  along 
the  z-axis  to  be  N  =  640.  The  spatial  step  size  in  both  the  x  and  z  directions  is  given 
by 

Ax  =  0.4/A.  (41) 

Then,  from  (40)  we  obtain  At  ~  1.0417  x  10-12  seconds. 

The  parameters  that  give  the  position  of  the  source  term  Js  defined  in  (9)  are 
x  =  0.025  and  z  =  0.1.  For  our  simulations,  the  frequency  of  our  sine  wave  is  3  GHz, 
and  hence  u  =  27r*3  x  109.  The  duration  tf  of  the  sine  pulse  is  0.667  nanoseconds  (ns), 
which  is  two  complete  periods. 

The  depth  of  the  PML  layers  is  0.05  meters.  This  measurement  balances 
computational  efficiency  with  the  minimization  of  numerical  noise  caused  by  the 
reflections  from  the  PML  interfaces.  The  spatial  discretization  within  the  PML  is  the 
same  as  that  within  the  computational  region. 

In  order  to  attempt  the  parameter  identification  problem,  we  must  collect  data 
within  the  computational  domain  at  a  pre-specihed  sampling  rate.  We  collect  data 
at  the  center  of  the  antenna  at  each  time  step.  Then,  our  data  consists  of  the  set 
{Ey(iAt,  0,  z,  q)}^1;  where  Ey  is  the  solution  of  Maxwell’s  equations  given  by  the  FDTD 
algorithm,  q  =  (a,  r,  es,  e^)  is  the  set  of  parameters  that  characterize  the  Debye  medium, 
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Figure  3.  Data  for  the  Debye  model  for  water.  The  vertical  axis  gives  the  magnitude 
of  the  electric  field.  The  horizontal  axis  gives  the  time  in  the  units  At  seconds. 


and  Nt  is  the  total  number  of  time  steps.  Notice  that  we  have  now  made  explicit  the 
dependence  of  solutions  of  (31)-(38)  on  the  parameter  set  q.  We  will  continue  to  use  this 
convention  in  the  sequel.  For  the  discretization  mentioned  above,  we  take  Nt  =  1800,  in 
which  case  the  data  includes  both  the  outgoing  energy  and  the  energy  that  is  reflected 
from  the  free-space/diclectric  interface. 

In  [5]  it  is  noted  that  the  polarization  behavior  of  water  is  reasonably  modelled  by 
(37),  (39),  and  (38)  with  parameter  values: 

a  =  1  x  10"5  mhos/meter, 

r  =  8.1  x  10~12  seconds,  (42) 

es  =80.1  relative  static  permittivity, 

Coo  =  5.5  relative  high  frequency  permittivity. 

With  this  information,  numerical  approximations  to  (31)-(38)  can  be  obtained.  The 
resulting  data  set  {Ey(iAt,  0,  z,  q)}^  is  plotted  in  Figure  3.  In  Figure  4,  we  see  the 
formation  of  the  Brillouin  precursors  [1,  2,  3]  within  the  dielectric.  This  suggests  that 
the  FDTD  solutions  are  accurately  capturing  the  dispersivity  of  the  Debye  medium. 
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Figure  4.  Precursor  Formation  Within  the  Debye  Medium.  We  plot  {Ey(t,  0,  z)},  for 
z  £  [0.2, 0.4]  at  four  different  times.  In  the  upper  left,  t  =  1500At.  In  the  upper  right, 
t  =  2500AI.  In  the  lower  left,  t  =  3500.  In  the  lower  right,  i  =  4500.  These  times 
correspond  to  1.56,  2.60,  3.65,  and  4.69  ns  respectively. 


3.  The  Inverse  Problem 

We  now  address  the  problem  of  parameter  identification  for  a  Debye  medium.  In  the 
previous  section,  the  problem  of  solving  Maxwell’s  equations  with  a  Debye  polarization 
model  was  addressed.  Implicit  in  this  problem  is  the  knowledge  of  the  “true”  parameter 
vector  q*,  given  by  (42).  In  the  inverse  problem,  on  the  other  hand,  we  have 
measurements  of  the  electric  field  at  the  point  (0,  z)  in  the  computational  domain  at 
uniform  intervals  of  time  tt.  Given  this  data,  our  goal  is  to  then  identify  the  “true” 
parameter  vector  q*. 

3.1.  The  Noise  Model  for  Data  Generation 

In  the  laboratory  setting,  noise  enters  into  the  problem  both  through  the  resistor,  which 
generates  the  electric  pulse,  and  through  the  antenna/receiver,  which  records  the  electric 
field  intensities  at  the  point  (0,  z)  at  discrete  time  steps.  The  noise  model  presented  in 
this  section  follows  [13]  closely. 

First,  noise  enters  the  signal  via  the  random  motion  of  electrons  in  the  resistor. 


Parameter  Identification  for  a  Dispersive  Dielectric 


13 


These  random  motions  produce  small,  random  voltage  fluctuations  at  the  resistor 
terminals.  Over  time  these  voltage  fluctuations  have  a  zero  average  value,  but  non 
zero  rms  (or  standard  deviation)  given,  at  time  tt,  by  Planck’s  black  body  radiation 
law, 


where 


Vi 


4  hfBR 

ehf/KT  _  i  ’ 


(43) 


h  =  6.546  x  1CT34  J-sec  is  Planck’s  constant, 
k  =  1.380  x  1CT23  J/°K  is  Boltzmann’s  constant, 
T  is  the  temperature  is  degrees  kelvin  (K), 

B  is  the  bandwidth  of  the  system  in  Hz, 

/  is  the  center  frequency  of  the  bandwidth  in  Hz, 
R  is  the  resistance. 


For  microwave  frequencies  we  have  hf  «  kT ,  and  hence, 
phf/KT  _  1  _  hf 

KT' 

Thus  equation  (43)  is  well  approximated  by 

n  =  VikTBR.  (44) 

This  is  the  approximation  commonly  used  in  microwave  work  [13].  The  noise  is  modelled 
as  independent  and  identically  distributed  (i.i.d.)  Gaussian  with  constant  variance  vf. 

If  we  replace  the  noisy  resistor  with  a  noiseless  resistor  together  with  a  voltage 
generator  [13],  then  we  have 

Vi  =  VkTB.  (45) 


Noise  also  enters  the  data  through  the  antenna  and  receiver.  The  noise  comes  from 
two  sources:  the  noise  due  to  the  antenna,  and  the  electronic  noise  in  the  receiver.  Both 
noises  are  assumed  to  be  white  Gaussian  with  standard  deviation  V{  =  \JkTB  [10,  13]. 

With  this  model,  the  noise  in  the  data  is  zero  mean,  additive  Gaussian,  with 
constant  variance  across  time.  That  is,  our  data  vector 


•p ydata  _ 


(  TT'data 

~  (A/,  1  ^ 


rpdata\ 

^y,Nt ) 


is  a  random  vector  with  normally  distributed  components 


(46) 


S*‘“~V(B„(«j,0,f,q-),E2)  for  i  =  1 . JVt,  (47) 

where  £2  is  calculated  using  the  above  approximations  for  the  variances  of  the  noise 
that  comes  from  the  resistor,  the  antenna,  and  the  receiver. 

To  generate  synthetic,  noisy  data,  we  use  the  FDTD  method,  with  the  specifics  of 
Section  2.4,  to  obtain  the  set  {Ey(ti,  0,  z,  q*)}^,  together  with  noise  model  (47).  We 
then  seek  the  maximum  likelihood  estimator  (MLE)  of  the  true  parameter  vector  q* 
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given  by  (42).  Given  the  noise  model  (47),  with  data  equal  to  n  realizations  {dj}f=1  of 
(46), (47),  the  MLE  qn  for  q*  is  given  by 

1  71 

qn  =  arg min  J (q),  where  J(q)  =  -  I iEy(°,  z>  q)  -  dj 111,  (48) 

qGC^  Z 

where  ~Eiy(f),z,  q)  =  (Ey(ti,  0,  z,  q), . . . ,  Ey(tNt,  0,  z,  q)),  and  Q  is  the  constraint  set  for 
the  parameter  values. 

3.2.  The  Optimization  Problem 

The  minimization  problem  (48)  is  solved  using  a  modified  Levenberg-Marquardt  method 
[11].  In  the  particular  implementation  that  we  employ,  an  ellipsoidal  trust  region 
approach  is  taken,  which  alleviates  difficulties  caused  by  large  differences  in  the 
magnitudes  of  the  parameters  in  the  vector  q.  In  particular,  each  variable  is  scaled 
to  be  of  roughly  the  same  magnitude,  namely  each  scaled  variable  is  approximately 
0(1)  near  q*.  Prior  knowledge  of  the  magnitudes  of  the  parameters  at  the  solution  can 
be  used  to  determine  the  scaling  factors.  The  trust-region  subproblem  is  solved  at  each 
outer  iteration  using  the  dogleg  method  (see  [11]).  In  addition,  we  enforce  the  physical 
constraints  cr>0,r>0,  es>l,  and  >  1,  which  determine  Q.  This  is  done  via 
the  projection  of  both  the  trust  region  and  of  the  search  directions  onto  the  feasible 
set  at  each  iteration.  The  convergence  properties  of  the  unconstrained  version  of  the 
algorithm  are  then  retained.  The  details  of  this  algorithm  will  be  presented  in  a  future 
paper. 

The  Jacobian  matrix,  which  is  used  both  for  the  computation  of  the  gradient  and  for 
the  construction  of  the  Gauss-Newton  approximation  to  the  Hessian,  is  approximated 
using  forward  differences.  Each  iteration  then  involves  at  most  five  function  evaluations: 
one  at  the  current  iterate,  and  one  for  each  of  the  four  parameters.  This  is  not  overly 
restrictive  for  our  problem,  but  for  problems  in  which  the  permittivity  and  conductivity 
are  allowed  to  be  spatially  varying,  this  is  no  longer  a  viable  approach.  Analytic 
computations  of  the  gradient  are  possible,  would  alleviate  this  difficulty,  and  require 
no  finite  difference  approximations.  We  hope  to  explore  the  use  of  analytic  gradients 
for  this  problem  in  future  work. 

It  is  important  to  note  that  a  solution  to  (48)  depends  upon  the  spatial 
discretization  Ax  given  by  (41)  (note  that  At  is  then  obtained  using  (40)),  which  depends 
only  on  the  number  of  grid  points  N  in  Zv  U  Zd  of  Figure  2.  (We  assume  that  the  length 
of  0.4  for  ZvUZd  is  fixed.)  Thus  in  (48),  J  and  any  approximate  solution  q  of  (48)  should 
carry  the  index  N.  Nonetheless,  in  the  sequel  we  will  suppress  this  N  dependence. 

3.3.  The  Inverse  Problem  for  a  Debye  Medium 

We  now  present  results  for  the  inverse  problem  of  identifying  the  dielectric  parameters 
in  the  Debye  model  for  water  given  noisy  data.  Using  the  parameters  given  in  (42) 
and  the  discretization  details  given  in  Section  2.4,  we  generate  data  with  varying  levels 
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Table  1.  Estimated  parameters  in  Debye  inverse  problem.  Test  1 


Test 

%  Noise 

a 

r 

<7 

^OO 

Residual 

True  values: 

1.0  x  10"5 

8.1  x  10-12 

80.1 

5.5 

1 

0.0 

0.0 

9.965  x  10-12 

80.75 

14.77 

6.7051  x  10"8 

2 

0.5 

0.0 

9.974  x  10”12 

80.20 

13.30 

1.2279  x  10“5 

3 

1.0 

0.0 

9.980  x  10"12 

80.18 

12.23 

4.8912  x  10"5 

4 

1.5 

1.952  x  10-1 

7.437  x  10”12 

79.07 

10.00 

1.1483  x  10~4 

5 

2.0 

0.0 

9.989  x  10-12 

80.01 

10.66 

1.9542  x  10"4 

6 

2.5 

3.671  x  10~2 

9.997  x  10~12 

79.57 

9.717 

3.1861  x  10"4 

7 

3.0 

8.715  x  10“2 

9.998  x  10”12 

79.82 

9.189 

4.5844  x  10~4 

8 

5.0 

0.0 

9.903  x  IQ”12 

85.42 

16.22 

1.2000  x  10“3 

of  noise  using  the  error  model  (47)  with  tl  =  iAt.  Hence,  our  sampling  rate  for  the 
antenna/receiver  is  assumed  to  be  At  «  0.001  ns.  Admittedly,  such  a  sampling  rate 
may  not  be  possible  in  many  laboratory  settings.  Several  numerical  experiments  were 
performed  in  which  data  was  generated  using  error  model  (47)  with  fj  =  i- 1  ■  At,  where 
I  is  an  integer  larger  than  1.  As  one  would  expect,  our  ability  to  identify  the  true 
parameter  vector  q*  in  the  inverse  problem  decreased  as  t  increased.  Once  real  data 
is  obtained,  with  a  particular  sampling  rate,  it  will  likely  be  necessary  to  decrease  the 
sampling  rate  when  solving  the  inverse  problem. 

We  make  eight  different  attempts  at  the  inverse  problem  with  a  different  noise  level 
for  each  attempt.  We  solve  (48)  using  the  approach  discussed  in  the  previous  section 
with  initial  guess  q0  =  (cr0,  r0,  eSi0,  Coo,o)  given  by  cr0  =  1.5  x  10~5,  r0  =  10.0  x  10~12, 
eSjo  =  73.1,  and  e^o  =  6.0.  (These  range  between  50%  relative  to  10%  relative  error 
from  the  true  values  and  are  typical  of  values  used  in  testing  algorithms  [4].).  We 
use  only  one  data  vector  di.  Then  n  —  1  in  (48).  The  results  are  given  in  Table  1. 
The  residual,  which  is  given  in  the  far  right  column,  is  J(q)/||di||2,  where  q  is  the 
approximate  solution  to  (48)  given  by  the  optimization  algorithm. 

A  partial  explanation  for  results  in  Table  1  can  be  obtained  from  an  analysis  of 
equation  (14),  which  we  rewrite  as 

*/  N  _  G  ,  J^T  ,  V  /„0N 

er(u>)  —  7“ — : - TV — : eoo  +  V •  (49) 

1+JUT  l+JOJT  JLU60 

Since,  in  our  simulations,  the  outgoing  radiation  is  given  by  two  cycles  of  a  sine 
wave,  in  the  continuum  space,  the  resulting  pulse  will  have  infinite  frequency  content. 
Nonetheless,  for  the  purpose  of  the  inverse  problem,  the  outgoing  and  reflected  radiation 
will  be  dominated  by  frequencies  near  the  carrier  frequency,  3  GHz.  This  is  seen  in  Figure 
5.  Thus  e*  will  be  dominated  by  frequencies  near  3  GHz,  and  so  we  restrict  our  analysis 
of  (49)  to  the  case  where  u  m  3  x  109.  Then,  since  lot*  ~  (T(1CT2), 


0{£ 


1  +  joJT 


r^j 


(50) 
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Figure  5.  Fourier  Transform  of  the  Data  in  Figure  3.  The  vertical  axis  gives 
amplitude.  The  horizontal  axis  gives  the  frequency  from  0.3  GHz  to  30  GHz. 


JUT*  * 

1  +  jUT*  600 


0(10 


— 2  * 


(51) 


Furthermore,  given  that  eo  ~  8  x  10  12,  we  have 

— - C?(10V).  (52) 

juto 

Equations  (50),  (51),  and  (52)  suggest  that  e*  will  be  most  sensitive  to  the  static 
permittivity  es.  More  specifically,  near  q*,  e*  will  is  more  sensitive  to  changes  in  es  than 
to  changes  in  r,  e^,  and  a  of  the  same  relative  magnitude.  Further  support  for  this 
claim  can  be  found  by  an  appeal  to  the  reflection  coefficient  T,  which  is  defined  by 


r 


1  —  y/e*r 


r  will  also  be  dominated  by  es.  We  therefore  expect  that  the  cost  function  J  will  be 
most  sensitive  to  es,  and  hence,  that  es  will  be  the  easiest  parameter  to  reconstruct. 
This  is  supported  by  the  results  in  Table  1. 

Furthermore,  equations  (50)  and  (51)  suggest  that  near  r*,  e*  is  relatively  insensitive 
to  “small”  changes  in  r,  since  if  r  ~  t*,  the  approximation  1  +  jojr  «  1  remains  valid. 
On  the  other  hand,  if  r  becomes  too  large,  it’s  effects  will  be  manifested  in  both  the  first 
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and  second  terms  of  (49).  In  fact,  we  found  that  if  r  is  fixed  and  is  large  in  comparison 
to  t*  ,  the  reconstructed  values  of  es  are  very  poor.  These  observations  explain  why  in 
Table  1  the  reconstructed  values  of  r  are  near  in  magnitude  to  t*  at  all  noise  levels. 

Since  e*  is  dominated  by  es,  we  would  not  expect  the  accurate  reconstruction  of  e^. 
Furthermore,  assuming 

0(t)  =  0(t *),  (53) 

re oo  «  T*e*ooi  (54) 

from  equation  (51),  we  have 

JUT  juT*  * 

— eoo  «  — -e  .  55 

1  +  JUT  1  +  JUT* 

Hence,  independent  identification  of  t*  and  e^,  will  be  very  difficult,  particularly  in  the 
presence  of  noise.  The  identification  of  their  product,  on  the  other  hand,  will  be  less 
difficult,  though  we  would  not  expect  to  reconstruct  r*^  as  accurately  as  we  do  e*. 
This  is  also  supported  by  the  results  of  Table  1. 

The  observations  of  the  previous  two  paragraphs  suggest  a  strategy  for  the  inverse 
problem.  Namely,  fix  the  value  of  e to  an  approximate  of  of  the  same  order  and 
optimize  over  the  other  parameters.  The  optimized  value  of  r  will  then  likely  satisfy 
(53),  (54),  and  (55),  in  which  case  the  dynamics  of  the  problem  will  be  adequately 
captured.  Before  testing  this  approach,  we  discuss  the  parameter  a. 

Equation  (52)  suggests  that  J  is  least  sensitive  to  the  parameter  a.  In  addition, 
since  the  change  in  cr  from  its  value  in  free  space  to  its  value  in  the  dielectric  is  very 
small,  we  would  expect  very  little  of  the  reflected  energy  to  be  due  to  the  change  in  a. 
The  results  of  Table  1  are  therefore  no  surprise.  Also  not  surprising  is  the  fact  that  the 
magnitude  of  o  remains  small  in  all  reconstructions.  This  suggests  the  strategy  of  fixing 
the  value  of  cr  to  a  value  reasonably  close  to  cr*  prior  to  optimization.  The  dynamics  of 
the  problem  should  then  be  adequately  captured. 

We  now  test  the  strategies  suggested  in  the  previous  two  paragraphs.  Namely,  we 
perform  the  inverse  problem  with  the  values  of  a  and  e fixed  at  cr  =  0.0  and  =  1.0. 
The  results  of  this  approach  are  presented  in  Table  2.  The  values  of  eSjo  and  To  are  the 
same  as  in  the  previous  attempt.  Comparing  these  results  with  those  of  Table  1,  we 
see  that  this  strategy  yields  very  little  change  in  the  values  of  the  residuals,  while  the 
reconstruction  of  es  remains  accurate.  In  addition,  the  solution  to  the  inverse  problem 
can  be  had  far  more  quickly  when  cr  and  e ^  are  fixed  at  the  outset.  Though  more 
iterations  are  needed  in  this  case,  gradient  computations  are  obtained  with  3/5ths  the 
cost.  Also,  once  cr  and  e ^  are  fixed,  constraints  are  no  longer  necessary. 

The  relative  stability  of  the  r  reconstruction  in  Table  2  suggests  yet  another 
strategy.  Namely,  fix  the  value  of  r  at  the  outset.  Results  for  this  approach  are  given 
in  Table  3  where  r  was  fixed  to  be  the  initial  guess  r0  from  the  previous  attempts. 
Surprisingly,  we  find  the  this  strategy  actually  yields  slightly  smaller  values  for  the 
residuals  than  were  obtained  in  Table  2.  In  addition,  the  optimization  problems 
converged  far  more  quickly.  In  fact,  at  all  noise  levels,  within  three  iterations  our 
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Table  2.  Estimated  parameters  in  Debye  inverse  problem.  Test  2. 


Test 

%  Noise 

a 

T 

es 

^OO 

Residual 

True  values: 

1.0  x  10~5 

8.1  x  10~12 

80.1 

5.5 

1 

0.0 

0.0 

9.818  x  lO’12 

80.82 

1.0 

5.8619  x  10~7 

2 

0.5 

0.0 

9.767  x  10”12 

80.42 

1.0 

1.3168  x  10~5 

3 

1.0 

0.0 

9.781  x  10"12 

80.94 

1.0 

5.2418  x  10"5 

4 

1.5 

0.0 

9.751  x  10”12 

82.09 

1.0 

1.1127  x  10"4 

5 

2.0 

0.0 

9.922  x  10~12 

82.69 

1.0 

1.9741  x  10“4 

6 

2.5 

0.0 

9.123  x  10~12 

79.82 

1.0 

3.1841  x  10~4 

7 

3.0 

0.0 

9.598  x  10”12 

79.44 

1.0 

4.5081  x  10"4 

8 

5.0 

0.0 

9.703  x  10”12 

79.51 

1.0 

1.2000  x  10"3 

Table  3.  Estimated  parameters  in  Debye  inverse 

problem 

.  Test  3. 

Test 

%  Noise 

o 

T 

es 

^OO 

Residual 

True  values: 

1.0  x  10~5 

8.1  x  lO’12 

80.1 

5.5 

1 

0.0 

0.0 

1.000  x  10“n 

80.75 

1.0 

6.8509  x  10~7 

2 

0.5 

0.0 

l.ooo  x  hr11 

81.24 

1.0 

1.2546  x  10“5 

3 

1.0 

0.0 

1.000  x  HT11 

80.03 

1.0 

5.0736  x  10"5 

4 

1.5 

0.0 

1.000  x  10”11 

80.26 

1.0 

1.1028  x  10"4 

5 

2.0 

0.0 

1.000  x  10~n 

82.46 

1.0 

1.9746  x  10“4 

6 

2.5 

0.0 

1.000  x  10~n 

80.11 

1.0 

3.1851  x  10-4 

7 

3.0 

0.0 

1.000  x  10”11 

81.74 

1.0 

4.4774  x  10~4 

8 

5.0 

0.0 

1.000  x  10”11 

78.72 

1.0 

1.2000  x  10"3 

optimization  algorithm  had  converged.  Using  the  previous  strategies,  many  more 
iterations  were  required. 

The  previous  approach  rested  heavily  on  the  fact  that  our  choice  of  To  was  near  in 
magnitude  to  r*.  An  obvious  question,  then,  is  what  can  be  done  if  one  cannot  obtain 
a  good  initial  guess  for  r*?  We  take  the  approach  of  fixing  the  parameters  cr  =  0, 
es  =  73.1,  and  =  1  and  optimizing,  with  no  noise,  with  respect  to  the  parameter  r. 
With  an  initial  guess  of  To  =  lO-4,  we  obtain  the  approximation  r  =  6.9987  x  1CU12.  We 
observe  how  close  this  is  to  r*.  We  then  fix  this  value  of  r  and  optimize  over  es.  The 
results  of  this  approach  are  given  in  Table  4  and  are  similar  to  those  found  in  Table  3. 

3-4-  Reconstructions  with  Inferior  Accuracy 

In  each  of  the  above  experiments,  the  data  was  generated  on  the  same  grid  that  was 
used  in  the  inverse  problem,  and  i.i.d.  Gaussian  noise  was  added.  In  this  section,  we 
attempt  the  inverse  problem  with  data  generated  instead  on  a  finer  grid.  No  additional 
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Table  4.  Estimated  parameters  in  Debye  inverse  problem.  Test  4. 


Test 

%  Noise 

a 

r 

^OO 

Residual 

True  values: 

1.0  x  10~5 

8.1  x  10-12 

80.1 

5.5 

1 

0.0 

0.0 

6.9987  x  10-12 

79.82 

1.0 

5.7768  x  10"8 

2 

0.5 

0.0 

6.9987  x  10-12 

79.67 

1.0 

1.2299  x  10“5 

3 

1.0 

0.0 

6.9987  x  10-12 

80.18 

1.0 

5.4274  x  10“5 

4 

1.5 

0.0 

6.9987  x  10"12 

80.02 

1.0 

1.1609  x  10"4 

5 

2.0 

0.0 

6.9987  x  10-12 

78.54 

1.0 

2.0476  x  10~4 

6 

2.5 

0.0 

6.9987  x  10-12 

78.64 

1.0 

3.1835  x  10~4 

7 

3.0 

0.0 

6.9987  x  10-12 

82.33 

1.0 

4.4185  x  10~4 

8 

5.0 

0.0 

6.9987  x  IQ”12 

78.67 

1.0 

1.3000  x  10~3 

Table  5.  Estimated  parameters  in  Debye  inverse  problem.  Test  5. 


Test 

a 

r 

^OO 

Residual 

True  values: 

1.0  x  10~5 

8.1  x  10”12 

80.1 

5.5 

1 

0.0 

6.9987  x  10-12 

72.22 

1.0 

3.8477  x  10"5 

2 

1.0  x  10~5 

8.1000  x  10-12 

72.55 

5.5 

3.6494  x  10~5 

noise  is  added.  More  specifically,  we  generate  our  data  on  the  grid  with  N  =  1280  in 
(41)  and  attempt  the  inverse  problem  on  the  grid  with  N  =  640.  In  our  first  experiment, 
we  attempt  the  inverse  problem  with  a ,  r,  and  fixed  to  the  values  found  in  Table 
4.  These  results  are  given  by  Test  1  in  Table  5.  The  reconstructed  value  of  es  is 
disappointingly  far  from  e*.  In  order  to  see  if  the  choice  of  fixed  parameters  is  the 
culprit,  we  solve  the  inverse  problem  with  the  a,  r,  and  fixed  at  the  true  values  a*, 
t*,  and  e^,  respectively.  There  is  very  little  difference  in  the  result,  which  is  given  by 
Test  2  in  Table  5. 

One  possible  explanation  for  these  results  is  that  the  FDTD  algorithm  has  not  quite 
converged  for  N  =  640,  even  though  plots  of  {Ey(ti,  0,  z)},  such  as  is  found  for  N  =  640 
in  Figure  3,  seem  to  suggest  that  convergence  has  occurred.  Unfortunately,  solving  the 
inverse  problem  on  the  grid  defined  by  N  —  1280  in  (41)  is  computationally  challenging. 

If  the  FDTD  algorithm  has  not  converged  for  N  =  640,  then  one  is  led  to  question 
the  results  and  the  conclusions  of  Section  3.3.  Fortunately,  numerical  experiments 
suggest  that  they  will  remain  unchanged  for  N  =  1280.  Also,  for  N  =  320  and  N  = 
160,  the  same  experiments  were  performed  with  very  similar  results  and  conclusions. 
Moreover,  the  results  and  conclusions  of  Section  3.3  are  very  similar  to  those  found  in 
[3]  for  the  analogous  one  dimensional  problem  using  FEM.  Hence,  we  conclude  that 
they  are  valid.  However,  the  above  comments  do  provide  a  caution  that  the  values  of 
parameter  estimates  may  depend  significantly  on  resolution  of  the  underlying  dynamical 
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systems. 

One  should  also  question  the  strategy  of  using  FDTD  when  solving  the  inverse 
problem.  The  FDTD  algorithm,  though  second  order  accurate  in  free  space,  looses 
that  accuracy  within  the  PMLs  and,  more  importantly,  within  the  dielectric  [12]. 
Furthermore,  FDTD  is  highly  vulnerable  to  error  at  the  free  space/dielectric  interface. 
Thus,  though  FDTD  is  unmatched  in  simplicity  of  implementation  for  this  problem,  it 
may  not  be  the  proper  forward  solver  to  use  if  the  goal  is  the  solution  of  the  inverse 
problem.  Clearly,  further  investigation  is  necessary  before  this  question  and  similar  ones 
are  resolved. 

4.  Estimating  Variability  in  Reconstructions 

In  the  previous  section,  we  obtain  reconstructions  of  the  true  parameter  vector  q*  by 
approximately  solving  (48)  when  n  —  1.  In  each  of  the  above  approaches,  the  cost 
function  J  in  depends  upon  a  data  vector  dx,  which  is  a  realization  of  the  random  vector 
(46), (47).  If  the  experiment  that  was  performed  to  obtain  di  is  performed  a  second 
time,  a  different  data  vector  d2  is  obtained,  which  is  also  a  realization  of  (46), (47).  If 
the  computations  in  the  previous  section  are  repeated  using  d2  in  place  of  d1;  the  cost 
function  J  will  be  different,  and  hence,  different  values  for  the  reconstructed  parameter 
vector  q  will  be  obtained.  This  illustrates  the  fact  that  the  reconstruction  q(d)  is  in  fact 
random  as  well.  A  natural  question  therefore  arises.  Namely,  what  sort  of  variability 
is  in  our  reconstruction?  Or,  put  in  another  way,  what  kind  of  confidence  can  we  have 
that  our  reconstruction  q  is  close  to  the  true  parameter  vector  q*? 

To  address  this  question,  we  will  make  use  of  large-sample  theory  from  statistics 
(see  Chapter  7  in  [9]).  In  addition,  we  will  restrict  our  attention  to  the  case  of  Table  4, 
where  the  parameters  cr,  r,  and  are  fixed.  Then  the  cost  function  J  in  (48)  depends 
only  on  es.  We  now  state  the  theorem  from  [9],  p.  469,  that  we  will  use. 

Theorem  3.1:  Let  Dl5,.  ,,Dn  be  independent  and  identically  distributed  random 
vectors  that  depend  on  a  parameter  es  and  have  probability  density  function  f€a. 
If  fes  satisfies  certain  regularity  conditions,  then  any  consistent  sequence  esn  = 
dsn(Dl,  •  •  •  ,  D„)  of  roots  of  the  likelihood  equation  satisfies 

v^(eSn-e:)^A7(0,l//(e:)),  (56) 

where  e*  is  the  “true”  parameter  value,  and  it  is  assumed  that 

log  (57) 

satisfies  0  <  J(e*)  <  oo. 

The  random  vectors  Di, . . . ,  Dn  in  Theorem  3.1  are  given  by 
Bl  =  Edyata  i  =  1, . . .  ,n, 
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where  E data  is  defined  by  (46),  (47),  and  are  therefore  independent  and  identically 
distributed. 

The  density  function  f€g  is  related  to  the  cost  function  J  in  (48)  via 

J (es)  =  — £2  ■  log  fes  +  c,  (58) 

where  c  is  a  constant.  Hence, 

!«)=£..  .  (59) 

In  order  to  use  Theorem  3.1,  we  must  assume  that  the  cost  function  J  has  a  bounded 
third  derivative.  This  amounts  to  the  assumption  that  the  mapping  from  es  to 
Ey(t,Q,z,€s)  has  a  bounded  third  derivative,  which  will  be  true  almost  everywhere. 
The  remaining  regularity  conditions  required  for  fes  by  Theorem  3.1  are  easily  satisfied. 

Numerical  experiments  suggest  that  J(es)  is  strictly  convex,  in  which  case  the  only 
local  minimum  of  J  is  the  global  minimum.  Then,  if  our  algorithm  creates  a  sequence  for 
which  the  derivative  of  J  converges  to  zero,  this  sequence  must  converge  to  the  maximum 
likelihood  estimator.  Given  that  this  occurs  for  each  n,  the  sequence  {Is,n}™=i  will  be 
consistent. 

From  (56)  we  know  that 

^  * 

£s,n  ^ 

We  define 

OLn  =  -^-J"{eStn).  (60) 

a  zn 

It  is  shown  in  [9]  that 

I{fs,n)  *  0, 

and  hence 

«n-/(e:).  (61) 

We  will  therefore  use  an  to  approximate  J(e*)  in  (57). 

In  Table  6  we  solve  problem  (48)  for  various  values  of  n  at  the  2%  noise  level.  We 
note  that  as  n  increases,  the  values  for  es<n  and  an  converge,  then  appear  to  begin  to 
diverge  at  n  =  64.  The  divergence  is  likely  due  to  the  fact  that  for  large  values  of  n  the 
finite  difference  approximations  to  the  gradient  and  Hessian  become  inaccurate  for  the 
function  J  in  (48). 

The  rapid  convergence  of  {an}  suggests  that  the  approximations  of  J(e*)  are 
reasonably  accurate  for  small  values  of  n,  and  hence,  that  the  values  obtained  for  the 
approximate  standard  deviation,  which  are  found  in  the  last  column  of  Table  6,  should 
be  reasonably  accurate  as  well.  This  is  supported  by  the  results  found  in  Table  6,  where 
we  see  that  for  each  n,  with  the  exception  of  n  —  1  and  n  =  64,  eS)n  is  within  one 
standard  deviation  of  the  true  value  e*.  For  n  —  1  and  n  =  64,  the  reconstruction  is 
within  two  standard  deviations  of  e*,  which  is  also  reasonable.  We  note  also  that,  as 
predicted  by  Theorem  3.1,  eSjn  converges,  more  or  less,  to  the  true  value  e*  =  80.1. 
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Table  6.  Estimates  for  the  Variability  in  es  for  the  Example  with  2%  Noise.  The 
values  of  a,  r,  and  eoo  are  fixed  at  the  values  found  in  Figure  4.  We  use  an  to  estimate 
/(e*)>  1  /nan  to  estimate  the  variance,  and  1  / y/nan  to  estimate  the  standard  deviation. 


n 

^ s,n 

Oin 

1  /nan 

1  /y/nOLn 

1 

82.07 

0.3906 

2.5604 

1.6001 

2 

80.36 

0.4146 

1.2060 

1.0982 

4 

79.84 

0.4211 

0.5937 

0.7705 

8 

80.38 

0.4138 

0.3021 

0.5496 

16 

80.11 

0.4177 

0.1496 

0.3868 

32 

80.01 

0.4194 

0.0745 

0.2729 

64 

79.73 

0.4233 

0.0369 

0.1921 

Finally,  we  note  that  for  n  —  1,  the  approximation  to  /(e*)  is  near  to  the 
converged  value.  Also,  the  value  of  eS)i  is  reasonably  close  to  the  true  value  e*. 
Hence,  the  corresponding  approximate  to  the  standard  deviation  may  give  a  reasonable 
approximation  of  the  true  variability  in  the  reconstruction  eS)i.  In  the  next  experiment 
we  test  this  hypothesis. 

In  Table  7  we  address  the  following  question:  If  we  have  collected  data  only  once,  i.e. 
n  —  1  in  Theorem  3.1,  which  is  often  the  case  in  applications,  can  the  asymptotic  result 
of  Theorem  3.1  provide  a  reasonable  estimate  for  the  variability  in  the  corresponding 
reconstruction  eS;i  of  e*?  We  perform  the  experiment  for  eight  noise  levels,  and  compute 
the  corresponding  maximum  likelihood  estimate  esp.  We  approximate  the  standard 
deviation  of  esp  as  in  Table  6. 

First,  we  note  that  as  the  noise  level  rises,  so  does  the  variability  in  the 
corresponding  reconstruction  es,i.  Furthermore,  for  each  noise  level,  the  reconstructions 
are  well  within  two  standard  deviations  of  the  true  value  e*.  (See  Table  4  for 
reconstructions  with  n  —  1  for  different  data  and  at  the  same  noise  levels.)  Most  of 
the  reconstructions,  in  fact,  are  near  to  or  are  within  one  standard  deviation  of  e*.  We 
conclude  therefore  that  for  this  particular  problem,  the  approximations  of  the  standard 
deviation  found  in  Table  7  provide  one  with  a  reasonable  notion  of  the  variability  in  the 
reconstructions  es,i- 

Finally,  from  the  above  information  we  can  obtain  confidence  intervals  for  our 
reconstruction  eSjfl.  For  this  we  use  the  Wald  test  (  p.  525  of  [9]),  which  says  that  under 
the  regularity  conditions  of  Theorem  3.1  we  have  that 


(is,n-<)Vnin^  W(0,1),  (62) 

where  In  is  a  consistent  estimator  of  /(e*).  Then 

U/3/2  *  „  .  U/3/2 

£s,n  I — —  ^  ^s,n  T  j — —  •  (63) 

ynln  y  nln 

are  confidence  intervals  for  e*  with  asymptotic  confidence  coefficient  1  —  /3,  and  up/ 2  is 
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Table  7.  Estimates  for  the  variability  in  es  for  n  =  1.  The  values  of  <r,  r,  and 
are  fixed  at  the  values  found  in  Figure  4.  We  use  an  to  estimate  7(e*).  Then  1  /nan 
estimates  the  variance  and  1  j yjnan  estimates  the  standard  deviation  in  (56). 


%  noise 

£s,n 

1  /nan 

1  /y/nan 

0.5 

80.37 

6.6239 

0.1510 

0.3885 

1.0 

80.93 

1.6242 

0.6157 

0.7847 

1.5 

81.50 

0.7085 

1.4115 

1.1881 

2.0 

82.07 

0.3906 

2.5604 

1.6001 

2.5 

82.65 

0.2447 

4.0867 

2.0216 

3.0 

83.23 

0.1664 

6.0092 

2.4514 

5.0 

85.63 

0.0550 

18.1739 

4.2631 

Table  8.  Estimates  for  the  Variability  in  es  for  the  Example  with  2%  Noise.  The 
values  of  a,  t,  and  are  fixed  at  the  values  found  in  Figure  4. 


n 

^ s,n 

1/  ^nan 

P  U/3/2 

yJnOLn 

e  +  U0/2 
ts’n  ^ 

1 

82.07 

1.6001 

77.95 

86.19 

2 

80.36 

1.0982 

77.53 

83.19 

4 

79.84 

0.7705 

77.86 

81.82 

8 

80.38 

0.5496 

78.96 

81.80 

16 

80.11 

0.3868 

79.11 

81.11 

32 

80.01 

0.2729 

79.31 

80.71 

64 

79.73 

0.1921 

79.23 

80.22 

the  upper  (5/2  point  of  the  standard  normal  distribution.  It  is  given  in  [9],  p.  526,  that 
an  in  (60)  is  a  consistent  estimator  of  /(e*),  and  hence  we  can  use  In  =  an  in  (63)  to 
obtain  confidence  interval  estimates. 

For  the  results  summarized  in  Table  6,  we  use  (63)  to  obtain  approximate  confidence 
intervals  for  (3  =  .99.  The  results  are  found  in  Table  8.  We  see  that  for  n  small, 
the  confidence  intervals  are  large,  but  as  n  becomes  large,  the  confidence  intervals 
become  small,  and  a  great  deal  more  certainty  can  be  inferred  about  the  value  of  the 
reconstruction.  We  note  that  (63)  is  an  asymptotic  result,  and  hence,  is  suspect  in  the 
case  of  small  n.  Nonetheless,  for  each  n  in  our  example,  the  true  value  e*  is  contained 
within  the  given  confidence  interval. 

For  the  results  summarized  in  Table  7,  we  again  use  (63)  to  obtain  approximate 
confidence  intervals  for  f3  =  .9.  The  results  are  found  in  Table  9.  The  true  value  e*  is 
contained  in  the  confidence  interval  for  each  noise  level.  Hence,  we  conclude  that  even 
for  n  —  1  the  confidence  intervals  given  by  the  Wald  test  (63)  give  a  reasonable  notion 
of  the  variability  in  the  reconstructions  of  e*.  Note  the  dramatic  increase  in  the  size  of 
the  confidence  intervals  as  the  noise  level  increases. 
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Table  9.  Estimates  for  the  variability  in  es  for  n  =  1.  The  values  of  <r,  r,  and 
are  fixed  at  the  values  found  in  Figure  4.  We  use  an  to  estimate  7(e*).  Then  l/nan 
estimates  the  variance  and  If  Jnan  estimates  the  standard  deviation  in  (56). 


%  noise 

£s,n 

1  /y/nan 

P  M/3/2 

S’n  1 JnOLn 

e  +  up/2 
ts’n  ^ 

0.5 

80.37 

0.3885 

79.73 

81.01 

1.0 

80.93 

0.7847 

79.64 

82.22 

1.5 

81.50 

1.1881 

79.55 

83.45 

2.0 

82.07 

1.6001 

79.44 

84.70 

2.5 

82.65 

2.0216 

79.32 

85.98 

3.0 

83.23 

2.4514 

79.20 

87.26 

5.0 

85.63 

4.2631 

78.62 

92.64 

5.  Conclusions 

We  have  presented  a  mathematical  model  for  two  dimensional,  time  domain,  TM  mode 
calculations  with  PML  absorbing  boundary  conditions,  in  the  presence  of  a  Debye 
medium.  The  resulting  set  of  PDEs  is  solved  using  the  FDTD  algorithm.  We  present 
numerical  results  from  forward  calculations  when  the  Debye  medium  parameters  are 
chosen  in  order  to  model  the  polarization  behavior  of  water. 

A  statistical  model  for  data  generation  is  presented,  and  we  solve  the  parameter 
identification  problem  by  obtaining  the  maximum  likelihood  estimate  corresponding  to 
the  data  and  noise  model.  For  us  this  involves  the  minimization  of  the  negative  log- 
likelihood  function,  which  takes  a  least  squares  form.  A  constrained,  ellipsoidal,  trust 
region  modification  of  the  Levenburg-Marquardt  algorithm  is  used  to  this  end. 

A  frequency  domain  analysis  is  presented  that  provides  an  explanation  of  results 
obtained  from  solutions  of  the  inverse  problem  for  data  generated  with  varying  levels 
of  noise.  Strategies  for  solving  the  inverse  problem  which  decrease  computational  cost 
and  are  based  on  these  conclusions  are  presented. 

Finally,  we  use  a  theorem  from  the  statistical  theory  of  large  samples  to  obtain 
estimates  of  the  variability  in  the  estimated  parameters.  Confidence  intervals  for  these 
estimated  parameters  are  also  presented. 
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